'''
Calculates nCk.
    n choose k = n! / [k! * (n-k)!]
'''

import operator as op
from functools import reduce
import math
import numpy as np


def nCk(n, k):
    k = min(k, n-k)
    numer = reduce(op.mul, range(n, n-k, -1), 1)
    denom = reduce(op.mul, range(1, k+1), 1)
    return numer // denom


# Or use:
def nCr(n,r):
    f = math.factorial
    return f(n) // f(r) // f(n-r)  # // avoids overflows


def kaoshimoni(n):  # n 为题库中题目数量
    p = 0
    for _ in range(1000):
        ss = list()
        for _ in range(20):
            s = list(range(n))
            s = np.random.permutation(s)
            ss.append(s[:20])
        flag = 0
        for k1 in range(20):
            s1 = ss[k1]  # 20 个选定的题目
            for k2 in range(k1+1, 20):  # 对于其他人
                s2 = ss[k2]  # 其他人选定的题目
                if len(set(s1).intersection(set(s2))) > 3:
                    flag = 1
                    break
            if flag == 1:
                p += 1
                break
    return 1 - p/1000  # 不出现重复的概率


import random

for x in range(1100, 5000):
    count = 0
    for q in range(0, 1000):
        l = []
        for t in range(0, 20):
            l.append(random.sample(range(0, x), 20))
        for a in range(0, 20):
            f = False
            for s in range(a + 1, 20):
                if len(set(l[a]) & set(l[s])) > 3:
                    f = True
                    break
            if f:
                count += 1
                break
    if count / 1000 < 0.05:
        print(count / 1000, x)
        exit()
    print(x, count / 1000)


if __name__ == "__main__":
    # 52 张牌，没有大小王，平均分给四个人，求至少一个人拿到至少两张A的概率。
    # 全部的分法：
    denom = nCk(52, 13) * nCk(39, 13) * nCk(26, 13) * nCk(13, 13) * math.factorial(4)
    # 取反情况：每个人都拿到了一张 A
    numer = nCk(48, 12) * nCk(36, 12) * nCk(24, 12) * nCk(12, 12) * math.factorial(4) * math.factorial(4)
    p_allGotA = numer / denom
    print("Probability: {:.3f}".format(1 - p_allGotA))  # 0.895
    
    
    # 20个人参加面试，每个人随机抽取20道题，请问试题库最少有多少道题？才能满
    # 足20人中任意两人相同的题目<=3的概率>=95%。
    N = 1310
    print(190 * nCk(20, 4) * (20/N)**4 - 0.05)
    print(kaoshimoni(1209))
    
    def getProblems(n):
        count = 0
        for _ in range(1000):  # 模拟 1000 词
            problems = []  # 20 * 20
            for _ in range(20):  # 选题目
                ps = random.sample(range(n), 20)
                problems.append(ps)
            
            flag = False
            for i in range(20):
                for j in range(i, 20):
                    if len(set(problems[i]).intersection(set(problems[j]))) > 3:
                        flag = True
                        break
                if flag:
                    count += 1
                    break
        
        return count / 1000
    
    for x in range(1100, 1300):
        p = getProblems(x)
        if p <= 0.05:
            print("Probability: {:.4f}".format(p))
